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ABSTRACT: 



Several finite differencing schemes are used to compute the natural 
frequencies and mode shapes of simply supported circular cylindrical 
shells in an attempt to determine the most accurate numerical model. 

Sets of difference equations are developed from the governing differential 
field equations and by minimizing the finite difference form of the 
Lagrangian energy function. Both finite differences and trigonometric 
expansions are used to model the circumferential behavior. Staggered 
or half-stations are used in addition to the conventional differencing 
schemes. The results indicate that the schemes using the trigonometric 
expansions are generally more accurate than those using finite differences 
for the circumferential derivatives. Furthermore, the conventional 
differencing scheme is shown to be as accurate as the half-station 
scheme when the field equation approach is used in conjunction with 
the trigonometric expansions . 



TABLE OF CONTENTS 



Page 

ABSTRACT 1 

TABLE OF CONTENTS 2 

INTRODUCTION 3 

DIFFERENTIAL EQUATIONS, NATURAL 

MODES, AND NATURAL FREQUENCIES 4 

NUMERICAL MODEIE 6 

Field Equation Approach 10 

Energy Approach 13 

COMPARISON 16 

CONCLUSIONS 22 

REFERENCES 24 

APPENDIX 25 

DISTRIBUTION LIST 32 



2 



VIBRATION ANALYSIS OF CYLINDRICAL 
SHELIS BY SEVERAL FINITE DIFFERENCE SCHEMES* 



by 

Robert E. Ball 
Associate Professor 

Naval Postgraduate School, Monterey, California 
INTRODUCTION 

There are several large-scale digital computer programs presently in 
operation that use the method of finite differences to predict the natural 
frequencies and mode shapes of general shell of revolution structures 
[ 1 , 2 , 3 ]. These programs use a variety of numerical schemes in the develop- 
ment of the governing finite difference equations . This paper uses several 
of these schemes, and some others, to compute the natural frequencies and 
modes of simply supported circular cylindrical shells. Donnell’s theory 
is employed, and the unknown dependent variables are the three displacements. 
Two fundamental approaches are considered in the development of the 
difference equations; one is an energy formulation in which the variation 
of the finite difference form of the Lagrangian energy function with respect 
to the displacement variables at each station is set equal to zero, the 
other uses finite difference approximations to the continuous derivatives 
in the governing field equations. Within these two approaches, both 
finite differences and trigonometric expansions are used to model the 
circumferential behavior. In addition, the three displacements are either 
defined at the same stations , as in the conventional finite difference 
method, or they are defined at different, or staggered, stations. The 
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latter technique is also known as the half-station method, and has been 
used in the static analysis of beam-columns [4] and cylindrical shell 
roofs [ 5 ], and in the vibration analysis of circular rings [6]. The 



station scheme yields solutions that are significantly more accurate than 
the solutions from the conventional, or whole -station, method. The goal 
of this paper is to determine the most accurate numerical model from 
among the several considered for predicting the natural frequencies and 
mode shapes of simply supported circular cylindrical shells. The 
applicability of that scheme to the general shell of revolution is also 
examined . 



Consider the thin, elastic, circular cylindrical shell of length L, 
thichness h, and radius a shown in Fig. 1. The positive directions of 
the axial coordinate x, the circumferential coordinate 0, the axial 
displacement U, the circumferential displacement V, and the normal dis- 
placement W are as indicated in Fig. 1. Donnell's partial differential 
equations for the vibration analysis of the cylinder can be obtained by 
minimizing the Lagrangian with respect to U, V, and W, and assuming the 
motion to be harmonic in time with a natural frequency 0. This leads to [ 7 ] 



results presented in references [4] - [6] indicate that the staggered 



DIFFERENTIAL EQUATIONS, NATURAL MODES, 
AND NATURAL FREQUENCIES 





/1+iA /1-u,. 



2 

0 ) v 




{—) u, §0 + (— ) v, §§ + v, ee + W, Q 



-uu 



,5 " v»e - w - 12 + 2w ’|§09 




where a subscript comma denotes differentiation, u = U/a, 




Young ' s 



4 




Figure 1. Geometry, Coordinates and Displacements of the 
Circular Cylindrical Shell 
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modulus, p is the mass density, and v is Poisson' s ratio. The boundary 
conditions for the simply supported cylinder are 



u,g = v = w = w, 5? =0 I = 0, i (2) 

where SL = L/a. 

The coordinate 9 can be removed from equations (l) by assuming 
u = u(§)cos n0, v = v(§) sin n0, and w = w(§) cos n0. This gives 



2 / 1— u \ /1+U\ 2 

u, § g - n (-g-) u + n (-£-) v,g + uw,g = - u> u 



/l+U-> ,1-Uv 

~ ^ ( 2 ^ "*" 2 



2 2 
nv-nw = -u>v 



► (3) 



^ p |i 2 

- uu,g - nv - w - - 2n‘ i w,-g + n w) = - u> w 



The natural modes of the cylinder that satisfy the simply supported 
boundary conditions (2) are 



u = u cos (mTT§/£) cos n© 
v = v sin (mnl/jfc) sin n© 
w = w sin (miT§/ & ) cos n0 

where u, v, and w are arbitrary constants. Substituting 
equations (l) leads to 






(4) 



equations (4) into 
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The coefficients a are given in Table I. 

pq. 

NUMERICAL MODELS 

Consider the finite difference network shown in Fig. 2a, where d is 
the nondimensional meridional increment and A is the circumferential angular 
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TABLE I Matrix Coefficients a , Equation (5) 
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P mc = (sinp m d)/d = (2sinp m d/2)/d 

(sin nA)/A T| = (2sin nA/2)/A 





(b) Model A (c) Model B 




(d) Model C (e) Model D (f) Model E 



Figure 2. Numerical Models 
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increment. A typical station is denoted by i, j, where i' and j are numbered 

such that 5- = id and 0. = jA. The conventional differencing scheme where 
i J 

u, v, and w are each defined at the same station is shown in Fig. 2b. A 
staggered station scheme where u, v, and w are each defined at different 
locations on the shell [5] is shown in Fig. 2c. The conventional scheme 
for the ordinary differential equations (3) is shown in Fig. 2d, and the 
half-station scheme corresponding to Fig. 2c is shown in Fig. 2e. The 
staggered scheme used in reference [2] with the energy approach is shown 
in Fig. 2f. These five schemes will be referred to henceforth as models 
A, B, C, D, and E respectively. 

The difference approximations to the partial derivatives for model A 
are of the form 



(I). . ■ 



(6a) 



i,J 



(^-f) = (f. , . - 2f. . + f .)/d 2 

i,j 1 1,0 1 »J 1_1 ’ J 

(ape). . = ^ f i+i,j+i " f i-i,j+i ' f i+i,j-i + f i-i,j-i ) ^ 2d ^ 

1 3 j 



(6b) 



(6c) 



0r) - - 4 *i + i,j + 6 f i,J - 4 f i-2,X «*> 



. • 

i,J 



Cr 2 ^ 2 ) [f i+l,j+l + f i-l,j+l + f i+l,j-l + f i-l,j-l' 2(f i+l,d + f i,J+l 

,4 



55 59 ' 



+ f . , . + f. . , ) + 4w. . J /d 
1-1 ,0 1,3-1' 1,3 ' 



(6e) 



The derivative approximations for model B are those given by equations 
(6b), (6d), (6e) and 
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(as), , (f i4, J " 



(6f) 



1,0 



(9^0).^ , ^ f i+l,j+i " f i,j+i " f i+l,J-| + 



(6g) 



2,0 



Corresponding approximations are used for the ordinary derivatives for 
models C and D. 



Field Equations Approach 

The finite difference form of equations (l) corresponding to Model A is 



(U, 



1+1,0 - 2u i,3 + Vi.jVa 2 + <¥> (“M+i - 2u i, j - Xj-iVS 



+ (2 r> (v i + l,3 + l - V i-1, j+1 - v i + l,J-l + Vl,i-l)/( to ) * ”(Vl,J • Vl >3 >/(2< 

i = 0,1,2,...! ; j = 1,2, ...J 



2 

= -O) u. 



1,0 



(7a) 



,1+u 



( 2 > (u i+l,j+l “ u i-l,j+l " U i+l,j-l + Vl^-l^^ + (v i+l,j " 2v i 



v^^.Vd 2 + (v. } . +1 - 2v i? . + v.^.^/A 2 + (w. 5 . +1 - w. 5 ._ 1 )/(2A) 



= v i,j i = 0,1,2,...! ; 0 = 1,2, . . . J 
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- <n >j+ i - - “i,J - TS K*2,J 

-‘‘"i+l.d + 6 “i,j - 4u i-i,o + w i-2,j )/a4 + 2 [“i+1,3+1 * * w i+l,M 

-i ^ 

+w. . , - 2(w. , • + W. . , + w. , . + w, ) + 4w. /(dA) 

i-1, j-1 v i+l,0 1,0+1 1-1,0 i,0-l 1>0J 

+ (w . „ - 4w. . , + 6w. . - 4w. . , + w. . J/A 4 1 
v i,0+2 1,0+1 i,0 i,0-l i,0- 2 J 



— _o> w. . i = 1,2, ...1-1 ; J — 1,2, .. .J 

1,0 



where I = X/d and J = 2tt/A. The stations i = 1+1 are ficticious stations 
located one increment off the ends of the shell. The finite difference 
form of the boundary conditions (2) for model A is 



u 



i+1,0 



- u. n . 

i-l,0 



V. . = w. . = w. , . 

1,0 1,0 1 + 1,0 



2w. . + w. , . = 0 

1,0 1 - 1,0 



i = 0,1; 



0 — 1 , 2 , . . J ( 7b ) 

Equations (7) constitute a set of 3(l-l)J + 12J homogeneous algebraic 

2 

equations that define an eigenvalue problem. The eigenvalues are -oj , and 
the eigenvectors are the natural modes. For this set of equations assume 
the modes to be 

u. . = u cos (mui/l) cos (2mT j/j) (8a) 

1.0 

v. . = v sin (mrri/l) sin (2nTTj/<J) (8b) 

1 . 0 

w. . = w sin ( m rri/l) cos (2nrro/j) (8c) 

1 5 3 

Note that these modes satisfy the boundary conditions (7b). Furthermore, 
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substituting equations (8) into equations (7a) eventually leads to equations 

(5) with coefficients a that are identical for every i and i defined in 

pq 

equations (7a). Hence, equations (8) are the correct modes for model A 

and would be the eigenvectors determined by a digital computer except for 

round-off error. The coefficients a for this model are given in Table I. 

P! 

A somewhat different set of difference equations is obtained when the 
half-station scheme of model B is used. The results are given in equations 
(Al) in the Appendix. The natural modes for this model are also given by 
equations (8), but with i replaced by i-H? in equation (8a) and j replaced 
by j+g- in equation (8b). As in the previous case, substituting equations 
(8) into equations (Al) leads to equations ( 5 ) for every i and j defined 
in equations (Al). The coefficients a are given in Table I. 

The same procedure can be followed using equations (3), where the 9 
dependence has been eliminated, and models C and D. The difference 
equations and boundary conditions for these two models are given in the 
Appendix as equations (A2) and (A3) respectively. The natural modes 
for both of these schemes are 



= u cos 


(mui/l) 


(9a) 


= v sin 


(mrri/l) 


(9b) 


= w sin 


(imri/l) 


(9c) 



except in model D, where i is replaced by i+g- in equation (9a). The co- 
efficients a for both models are given in Table I. 

pq 

Equations (A2) are analogous to the set of difference equations 
developed by the computer program SALCRS [3]- The two sets of equations 
are not identical because Sanders ' shell theory is used in SALORS', and 
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because the governing equations are four second order differential 
equations in terms of u, v, w, and m^, the nondimens ional meridional 
bending moment. However, it can be shown that by eliminating the discrete 
variable m^ from the difference equations, the four equation formulation 
used in SALORS reduces to the identical three equation formulation given 
by equations (A2), when Donnell's theory is used. 

Energy Approach 

Sets of difference equations can also be developed by equating to 

zero the variation of the finite difference form of the Lagrangian with 

respect to the discrete displacement variables and the Lagrangian 

multipliers. The Lagrangian L can be expressed in the form 

8 

L=n-T+Z X. 6 . (10) 

k=l K * 

where II is the shell strain energy, T is the shell kinetic energy, is 
a Lagrangian multiplier, and 6^ is a constraint condition. Chuang and 
Veletsos [ 5 ] have shown that when model A is used the difference equations 
developed from the energy approach are not the same as equations (7a) , 
which were developed from the field equation approach. However, when 
model B is used, the two approaches yield identical equations. This 
section will consider the development of the difference equations by the 
energy approach using the models C, D, and E where the derivatives with 
respect to 0 have been replaced with multiples of n. 

"fch 

The shell strain energy in the Donnell theory in the n mode can be 
shown to be 

n - / < n A + v (iia) 

2(l-v ) o 
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where 



2 2 
n A = + 2v ( nv_w ) u >g + ( nv_w ) 

L a 2 r, x2 0 2 4 2, 

+ ^2 L ^ w ’ 5 |' ) " 2vn ^’§1 + n w ] 



Ilg = (^) (v, ? - nu) 2 + (l-v) n 2 (w, ? ) 2 



The shell kinetic energy in the n 



th 



mode is taken as 



1 = l (t a + V d ? 

2(l-v ) o 

where 

T A = ( v , T ) 2 + ( w , T ) 2 

T B ' <“•/ 



(lib) 



and 



t — t J E/pd-v 2 ) /a where t is time. 



Consider first model C. The integrals in equation (11) can be 
approximated by 

I 

2 ~ I 



n = 



rra Eh 



2 (l-v ) i=0 

I 



(I1 A + Vi d i 



T = 



2 

rra Eh 



X (t. + 



2(l-v ) i=0 



< T A + Vi d i 



(12a) 



(12b) 



where (II. + 

expressed at the i^* 1 station respectively, and d^ = d except at i = 0 and 
i = I where d^ = d/2. The constraint conditions for the simply supported 
shell are v^ = w^ = v^_ = w = 0. Carrying out the variation of the Lagrangia 



11^)^ and (T^ + T^)^ represent the strain and kinetic energy 



l4 



with respect to the variables u , v , w , u , v v w_ , 

X 1? \ 2 , and I4 and assuming harmonic motion leads to equations (A4) 

in the Appendix. A comparison of equations (A4) with equations (A2) 

reveals several differences. In the recursion equations, the differences 

appear in the approximation to the second derivative, and the energy 

2 2 2 

approach introduces the additional term (a vn d /24) (w^^ - 4w^ + ^ + 6w^ 

- 4w i ^ + w^ _ 2 )/d 4 in equation (AUc ) . Perhaps a more significant difference 
is the fact that the five special equations (A4e) - (A4i) are dissimilar to 
the equations obtained by applying the boundary conditions (A4d) to equations 
(A4a) - (A4b) at i=0 and i=l. The consequence of this is that the modes 
given by equations (9) are not the eigenvectors of equations (A4). The 
fourth boundary condition in equation (A4d) also eliminates equation (9c) 
as the mode for w. 

For model D, II and T are evaluated at i, as in the previous scheme, 

A n 

and Tg are evaluated at i+|r. Hence 



but Ilg 



n = 



T = 



rra^Eh 

2(l-v 2 ) 



2 

Tta Eh 
2(l-v 2 ) 






1-1 




+ 



I 



i=0 




(13a) 



(13b) 



The constraints are v Q = w Q = v^ = w^ = 0. Performing the variation of the 
Lagrangian and assuming harmonic motion yield a set of difference equations 
identical to equations (Al). 

For model E ? equations (12) are used in conjunction with equations (6f) 

for (u,^).^ and ( v, ^i’ equali 011 (6a) for (w,^)^; equation (6b) for 

(w, ) ; u. = (u ._i + u i)/2, and v = (v.^ + v. i.)/2. The difference 

1 1 1 "T 1 ~3 1 1- 2 
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equations for this scheme are given by equations (A5) in the Appendix. 

Note that, as is the situation with model C, the three special equations 
(A5f) - (A5h) are dissimular to the equations obtained by applying the 
boundary conditions (A5d) to equations (A5a) - (A5c). Hence, the modes 
given by equations ( 9 ) are not the eigenvectors of model E. 

COMPARISON 

The natural frequencies from the differential equations and models 
A-D have been computed for the four shells a = .01, & = 2; a = .01, i - 10; 
a = .1, i = 2; and a - .1, l = 10; for several values of m and n. Three 
axial increment sizes d = l/20, J0./4O , and £/l00; and two circumferential 
increment sizes A = 2rr/l0 and 2tt/20 were used. Poisson's ratio was .3* 

Since there are three degrees of freedom per mode, there are three 
eigenvalues for each value of m and n. The results for the lowest of the 
three natural frequencies for the two shells a = .01, i = 10, and a = .1, 

1 = 2 are given in Tables IIA and IIB respectively. Similar results were 
obtained for the other two shells. The natural frequencies from three 
dimensional elasticity theory are also given in Tables IIA and IIB [8]. 
Rather than present the natural frequencies of the numerical models, the 
quantity 100(l-<u/u) ) is presented in Tables IIA and IIB where ou e is the 
natural frequency associated with the differential equations (l). Thus, 
the smaller the number, the more accurate the natural frequency. 

Two other measures of the accuracy of the numerical models have been 
computed. These measures are 
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TABLE IIA Natural Frequencies 
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Note: The value 0 is used when 100(1-10, 



TABLE IIB Natural Frequencies 
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results were obtained for the other two shells . 
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CONCLUSIONS 



An examination of the results given in Tables IIA, IIB, IIIA, and 
IIIB reveals the following: 

(1) Model A yields extremely inaccurate natural frequencies for 
all m and n > 0 for the longer, thinner cylinder. The results 
are particular poor for m = 1. The frequencies are more accurate 
for the shorter, thicker cylinder. 

(2) Model B yields considerably more accurate natural frequencies 
than model A for the longer, thinner shell, and about equally 
accurate frequencies for the shorter thicker cylinder, except 

for m = 1 where model B is the more accurate. In general, model 
B is the more accurate model for all four shells . 

( 3 ) Models C and D yield essentially the same natural frequencies 
for all modes for all four shells. In some instances, model C 

is slightly more accurate, and in others, model D is the more 
accurate. In general, models C and D yield more accurate 
frequencies than models A and B. 

Hence, model B is preferred over model A, and models C and D are 
preferred over models A and B, provided the field equation approach 
is used to develop the difference equations for model C. The conclusion 
that model B is more accurate than model A is in agreement with the 
results presented in reference ( 5 ). The conclusion that models C 
and D are of equal accuracy is somewhat surprising, in view of the fact 
that model D yields identical difference equations from both the field 
equation and the energy approach. 
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The staggered station schemes of models B and D appear to be 
ideally suited to Donnell's theory with the boundary conditions given 
by equation (2) because only even or odd ordered derivatives of the 
dependent variables appear in any one field or boundary equation. 
General shells of revolution with meridional curvature do not have the 
consistant even and odd ordering of derivatives. Hence, the staggered 
station schemes do not appear to be more suitable for the general 
shell of revolution than for the cylinder. On the other hand, the 
meridional curvature may have some influence. 
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